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We introduce a combined fully three-dimensional macroscopic/microscopic transport approach 
employing relativistic 3D-hydrodynamics for the early, dense, deconfined stage of the reaction and a 
microscopic non-equilibrium model for the later hadronic stage where the equilibrium assumptions 
are not valid anymore. Within this approach we study the dynamics of hot, bulk QCD matter, 
which is being created in ultra-relativistic heavy ion collisions at RHIC. Our approach is capable of 
self-consistently calculating the freezeout of the hadronic system, while accounting for the collective 
flow on the hadronization hypersurface generated by the QGP expansion. In particular, we perform 
a detailed analysis of the reaction dynamics, hadronic freezeout, and transverse flow. 
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I. INTRODUCTION 

A major goal of colliding heavy-ions at relativistic en- 
ergies is to heat up a small region of space-time to tem- 
peratures as high as are thought to have occurred dur- 
ing the early evolution of the universe, a few microsec- 
onds after the big bang In ultra-relativistic heavy- 
ion collisions, such as are currently being explored at 
the Relativistic Heavy-Ion Collider (RHIC), the four- 
volume of hot and dense matter, with temperatures above 
- 150 MeV, is on the order of ~ (10 fm)'*. The state of 
strongly interacting matter at such high temperatures (or 
density of quanta) is usually called quark-gluon plasma 
(QGP) 5j. 



The first five years of RHIC operations at ^Jsnn = 
130 GeV and ^Jsnn = 200 GeV have yielded a vast 
amount of interesting and sometimes surprising results 
1^ , many of which have not yet been fully evaluated 
or understood by theory. There exists mounting evidence 
that RHIC has created a hot and dense state of decon- 
fined QCD matter with properties similar to that of an 
ideal fluid 0, Q ~ this state of matter has been termed 
the strongly interacting Quark- Gluon- Plasma (sQGP). 

Heavy-Ion collisions at RHIC involve several distinct 
reaction stages, starting from the two initial ground 
states of the colliding nuclei, followed by the high den- 
sity phase in which a sQGP is formed up to the final 
freeze-out of hadrons. 

The central problem in the study of the sQGP is that 
the deconfined quanta of a sQGP are not directly ob- 
servable due to the fundamental confining property of 
the physical QCD vacuum. If we could see free quarks 



and gluons (as in ordinary plasmas) it would be trivial 
to verify the QCD prediction of the QGP state. How- 
ever, nature chooses to hide those constituents within 
the confines of color neutral composite many body sys- 
tems - hadrons. One of the main tasks in relativistic 
heavy-ion research is to find clear and unambiguous con- 
nections between the transient (partonic) plasma state 
and the observable hadronic final state (for a review on 
QGP signatures, please see i). 

One particular approach to this problem is the appli- 
cation of transport theory. Transport theory ultimately 
aims at casting the entire time evolution of the heavy-ion 
reaction - from its initial state to freeze-out - into one 
consistent framework. By tuning the physical param- 
eters of the transport calculation to data one can then 
infer from these the properties of the hot and dense QCD 
matter of the sQGP and compare these to the predictions 
made by Lattice Gauge Theory (LGT). 



SPECIFIC MODEL FOR HIGH-ENERGY 
HEAVY-ION COLLISIONS 



II. 



Relativistic Fluid Dynamics (RFD, see e.g. j9l.lialllj) 

is ideally suited for the QGP and hydrodynamic expan- 
sion reaction phase, but breaks down in the later, dilute, 
stages of the reaction when the mean free paths of the 
hadrons become large and flavor degrees of freedom are 
important. The most important advantage of RFD is 
that it directly incorporates an equation of state as in- 
put and thus is so far the only dynamical model in which 
a phase transition can explicitly be incorporated. In the 
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ideal fluid approximation (i.e. neglecting ofF-equilibrium 
effects) - and once an initial condition has been specified 
- the EoS is the only input to the equations of motion 
and relates directly to properties of the matter under 
consideration. The hydrodynamic description has been 
very successful 0, EE Q in describing the collective 
behavior of soft particle production at RHIC. 

Conventional RFD calculations need to assume a 
freezeout temperature at which the hydrodynamic evolu- 
tion is terminated and a transition from the zero mean- 
free-path approximation of a hydrodynamic approach to 
the infinite mean-free-path of free streaming particles 
takes place. The freezeout temperature usually is a free 
parameter which (within reasonable constraints) can be 
fitted to measured hadron spectra. 

The reach of RFD can be extended and the problem 
of having to terminate the calculation at a fixed freeze- 
out temperature can be overcome by combining the RFD 
calculation with a microscopic hadronic cascade model - 
this kind of hybrid approach (dubbed hydro plus micro) 
was pioneered i n 1151 a nd has been now also taken up 
by other groups [id ■ Its key advantages are that the 
freezeout now occurs naturally as a result of the micro- 
scopic evolution and that flavor degrees of freedom are 
treated explicitly through the hadronic cross sections of 
the microscopic transport. Due to the Boltzmann equa- 
tion being the basis of the microscopic calculation in 
the hadronic phase, viscous corrections for the hadronic 
phase are by default included in the approach. 

Here, we combine the hydrodynamic approach with 
the microscopic Ultra-relativistic Quantum-Molecular- 
Dynamics (UrQMD) model in order to provide an 
improved description of the later, purely hadronic stages 
of the reaction. Such hybrid macro/micro transport cal- 
culations are to date the most successful approaches for 
describing the soft physics at RHIC. The biggest advan- 
tage of the RFD part of the calculation is that it directly 
incorporates an equation of state as input - one of its 
largest limitations is that it requires thermalized initial 
conditions and one is not able to do an ab-initio calcula- 
tion. 



A. Hydrodynamics 

In the present paper we shall use a fully three- 
dimensional hydrodynamic model for the description 
of RHIC physics, especially focusing on Au -|- Au colli- 
sions at RHIC energies {^/snn = 200 GeV per nucleon- 



nucleon pair). Our original code for solving the hydrody- 
namic equations, which based on Cartesian coordinates 
01, has been modified to the description on the coor- 
dinate by longitudinal proper time r — y/T^~~z^ and 
1] = ^ ln[(^ + z)/(t — z)], in order to optimize the hydro- 
dynamic expressions for ultra-relativistic heavy-ion colli- 
sions. 

In hydrodynamic models, the starting point is the rel- 
ativistic hydrodynamic equation 



(1) 



where T'^'^ is the energy momentum tensor which is given 
by 



(2) 



Here e, p, U and g^" are energy density, pressure, four 
velocity and metric tensor, respectively. We solve the 
relativistic hydrodynamic equation Eq. ^ numerically 
with baryon number ub conservation 



(3) 



In order to rewrite the relativistic hydrodynamic equa- 
tion Eq. 1^ in the coordinate (t, x,y,ri), we introduce 
the following variables [i^ , 



= cosh Yl / cosh{YL ~ 7j), 
Vy = Vy cosh Yl I cosh(yL — ??) , 
Vr^ = tanh(Yz:, - 77), 



(4) 



where = iln((l-hw^)/(l-w^)), 77 = l/21n((i + z)/(t- 
z)). Equation in the explicit way, is rewritten in the 
Appendix. 

In order to solve the relativistic hydrodynamic equa- 
tions, we adopt Lagrangian hydrodynamics. In La- 
grangian hydrodynamics, the coordinates of the individ- 
ual cells do not remain fixed, but move along the fiux of 
the fluid. In the absence of turbulence during the expan- 
sion, Lagrangian hydrodynamics has several advantages 
over the conventional Eulerian approach: 

• computational expediency: a fixed number of cells 
can be utilized through the entire calculation. A 
Lagrangian hydrodynamic code can thus easily be 
employed even at ultra-high energy collisions such 
as at the Large Hadron Collider (LHC) where a 
large difference of scale exists between the initial 
state and the final state due to the large gamma 
factor and rapid expansion of the QCD matter. 
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• analysis efficiency: the adiabatic path of each vol- 
ume element of fluid can be traced in the phase 
diagram, making it possible to directly discuss the 
effects of the phase transition on physical observ- 
ables 0. 

Our algorithm for solving the relativistic hydrody- 
namic equation in 3D is based on the conservation laws 
for entropy and baryon number. Further details concern- 
ing the numerical method can be found in ref. (igf. 

B. Equation of State 

To solve the relativistic hydrodynamic equation, an 
equation of state (EoS) needs to be specified. The inclu- 
sion of an equation of state as input is one of the biggest 
advantages of RFD, which is so far the only dynamical 
model in which a phase transition can explicitly be incor- 
porated. In the ideal fluid approximation (i.e. neglecting 
off-equilibrium effects), the EoS is the only input to the 
equations of motion and relates directly to properties of 
the matter under consideration. Once the EoS has been 
fixed (e.g. through a lattice-QCD calculation) a compar- 
ison to data can be used to extract information on the 
initial conditions of the hydrodynamic calculation 21]. 

Lattice-QCD (IQCD) offers the only rigorous approach 
for determining the EoS of QCD matter. Calculations 
at vanishing baryon chemical potential suggest that for 
physical values of the quark masses (two light {u, d)- 
quarks and a heavier s-quark) the deconfinement tran- 
sition is a rapid cross-over rather than a first order phase 
transition with singularities in the bulk thermodynamic 
observables The critical temperature at = for 
the rapid cross over in the (2+1) flavor case was recently 
predicted to be = 172 ± 11 MeV 0. 

However, many QCD motivated calculations for low 
temperatures and high baryon densities exhibit a strong 
first order phase transition (with a phase coexistence re- 
gion) These two limiting cases suggest that there 
exists a critical point (second order phase transition) 0| 
at the end point of a line of first order phase transi- 
tions. Recently, the exploration of the phase diagram 
for large temperatures and small, but non- vanishing, 
values of the baryon chemical potential became possi- 
ble through the application of novel techniques, such as 
Ferrenbe rg-S wendsen re- weighting Taylor series ex- 
pansions 123, l29l o r simulations with an imaginary chem- 
ical potential |3d l3l] | . Although these techniques have 
allowed for considerable improvements, the location of 



the critical point in the T - fiB plane still has large 
theoretical uncertainties, due to the sensitivity to the 
quark masses and the lattice sizes used in the calcula- 
tions. The predicted value of / j*^ varies between 

^endpouU^j.^ -1-3 liiEiH, i.e. ^''^dpo^nt ^^^q _ 

420 MeV. 

For the calculation presented in this work, we use a 
simple equation of state with the first order phase tran- 
sition, which will allow us to compare our results to pre- 
vious hydrodynamic and hybrid calculations employing 
(1-1-1) dimensional 0| and (2-1-1) dimensional hy- 
drodynamic models. 

Above the critical temperature (Tc = 160 MeV at /i = 
MeV) , the thermodynamical quantities are assumed to 
be determined by a QGP which is dominated by massless 
u, d, s quarks and gluons. The pressure in QGP phase is 
given by 



PQ 



(32 + 21AVh!^4 



180 

2 



(5) 



where Nj is 3 and B is the Bag constant |35l I36| . For 
the hadronic phase we use a hadr on g as equation of state 
with excluded volume correction |37|. Here, the pressure 
for fermions is given by 

pT\T,{f^^}) = ^pj^^^H7^,A*.-Vbp"'^'Hr,{A^J)) 

i 

= (6) 

i 

where is the pressure of ideal hadron gas and Vq is 

excluded volume of hadrons whose radii are fixed to 0.7 
fm. In the low-temperature region the well-established 
(strange and non-strange) hadrons up to masses of 
2 GeV are included in the EoS (see tables Tab. UandllTI 
for a detailed listing). Although heavy states are rare in 
thermodynamical equilibrium, they have a larger entropy 
per particle than light states, and therefore have consid- 
erable impact on the evolution. In particular, hadroniza- 
tion is significantly faster as compared to the case where 
the hadron gas consists of light mesons only (see the dis- 
cussion in SElEilillilS) 

The hadronic states used in the EoS of our hydrody- 
namic calculation are identical to those used in the mi- 
croscopic model employed for the hadronic stage of the 
reaction (UrQMD, see section Hi E|l . This is necessary to 
ensure consistency between the properties of the hadron 
gas in the hydrodynamic as well as in the microscopic 
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TABLE I: Baryons and baryon-resonances treated in the 
model. The corresponding antibaryon states are included as 
well. 
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TABLE IL Mesons and meson-resonances, sorted with re- 
spect to spin and parity, treated in the model. 



picture and allows for a smooth transition from one de- 
scription to the other. UrQMD additionally assumes a 
continuum of color-singlet states called "strings" above 
the m ~ 2 GeV threshold to model 2 ^ n processes and 
inelastic processes at high CM-energy. For example, the 
annihilation of an p on an is described as excitation 
of two strings with the same quantum numbers as the 
incoming hadrons, respectively, which are subsequently 
mapped on known hadronic states according to a frag- 
mentation scheme. Since we shall be interested in the dy- 
namics of the f2-baryons emerging from the hadronization 
of the QGP, it is unavoidable to treat string-formation. 
The fact that string degrees of freedom are not taken 
into account in the EoS © does not represent a problem 
in our case since we focus on rapidly expanding systems 
where those degrees of freedom can not equilibrate . 

The phase coexistence region is constructed employing 
Gibbs' conditions of phase equilibrium. The bag param- 
eter of B = 385 MeV/fm'^ is chosen to yield the critical 



temperature Tc = 160 MeV at = 0. In the coexistence 
region of QGP phase and hadron phase (i.e. the mixed 
phase) we introduce the fraction of the volume of the 
QGP phase, X{Xfj_) (0 < A < 1) and parameterize energy 
density and baryon number density as 

eM(A,T*) = XeQ{T*) + {l-X}eHiT*), 
nBM{X,T*) = XnBQ{T*) + {l-X)nBH{T*), (7) 

where T* — T*{^i) is the value of temperature on the 
phase boundary which is determined by Gibbs' condi- 
tions of phase equilibrium. 




160 
T [MeV] 



FIG. 1: Entropy density as a function of temperature at /i = 
0. The critical temperature is 160 MeV. 

In a forthcoming publication, we will discuss the EoS 
dependence of physical observables utilizing a realis- 
tic IQCD based equation of state, i.e. an EoS with a 
crossover phase transition at high T and low /x, includ- 
ing the QCD critical point 44]. We are currently in 
the process of constructing such an EoS, which will cor- 
rectly capture critical phenomena around the QCD crit- 
ical point 01 • A simple parameterization of the EoS 
around the QCD critical point as presented in 4,5j is un- 
fortunately insufficient for a description of these phenom- 
ena, even though it already provides a marked improve- 
ment over currently used equations of state. 



C. Initial Conditions 

The initial conditions for the hydrodynamic calculation 
need to be determined either by adjusting an appropriate 
parametrization to data or by utilizing other microscopic 
transport model predictions for the early non-equilibrium 
phase of the heavy-ion reaction. 

Numerous studies exist for finding the appropriate ini- 
tial conditions for hydrodynamic models |39|, |4y, |43, |43, 
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01 - usually such an initial condition is given by the 
parameterization of the spatial distribution of the en- 
ergy or entropy density and baryon number density at 
an initial time tq. A comparison between final particle 
distributions calculated by the hydrodynamic model and 
experimental data can then be utilized to fix the values 
of the parameters for the initial conditions. However, 
this ansatz is problematic if no experimental data exist 
to tunc the initial conditions. Furthermore, one looses 
predictive and analytic power by treating the quantities 
governing the initial conditions as free parameters. Re- 
cently there have been several attempts to determine 
a set of initial conditions not from parameterizations 
and comparison to data, but via a calculation using the 
color glass condensate (CGC) model for the initial state 
I47I ksf as well as an approach combining perturbative 
QCD and the saturation picture 0|. A study of el- 
liptic flow by Hirano et al. has shown that additional 
dissipation during the early QGP stage is required if an 
initial condition based on the CGC |43| is used. Another 
interesting fact which is found in hydrodynamic analy- 
ses at RHIC is that thermalization is achieved on very 
short timescales after the full overlap of the colliding nu- 
clei: none of the hydrodynamic calculations which have 
successfully addressed RHIC data at the top energy of 
y/SNN = 200 GeV/nucleon have initial times later than 
To ^ 0.6 fm js^j^l- The physics processes leading to 
such a rapid thermalization have yet to be unambigu- 
ously identified js^ - note that hydrodynamics itself can- 
not address the question of thermalization, since it relies 
on the assumption of matter being in local thermal equi- 
librium. 

For our calculation we use a simple initial condi- 
tion which is parameterized based on a combination of 
wounded nucleon and binary collision scaling 0, [sj |3 . 
Similar parameterizations have been used in various hy- 
drodynamic models which have been successful in ex- 
plainin g nu merous experimental observations at RHIC 
haMml- We have chosen this common initial condi- 
tion for our investigation in order to describe the general 
features of our model and provide a base-line comparison 
to previous calculations under similar assumptions ~ in a 
subsequent publication we shall investigate the sensitiv- 
ity of our results to particular variations and assumptions 
regarding the choice of the initial conditions. 

We factorize the energy density and baryon number 
density distributions into longitudinal direction {H(r])) 



and the transverse plane {W{x, y; &)), which are given by 



nBix,y,ri) = nBmiixW{x,y;b)H{r]), 



(8) 



where fmax and n^max sue parameters which are maxi- 
mum values of energy density and baryon number den- 
sity. The longitudinal distribution is parameterized by 



H{rj) = exp [-(hi - r;o)/2a2] e{\Tj\ ~ ,70), 



(9) 



where parameters 770 and ct,, are determined by com- 
parison with experimental data of single particle dis- 
tributions. The function W{x, y; b) on the transverse 
plane is determined by the superposition of wounded 
nucleon scaling which is characteristic of "soft" particle 
production processes and binary collision scaling which 
is characteristic of "hard" particle production processes 
[fi^f . This function is normahzed by W^(0,0;0). In the 
wounded nucleon scaling, the density of wounded nucle- 
ons is given by 



(10) 



where bA = s + b ■ {bs = s — b ■ e^), cr is the total 
nucleon-nucleon cross section at Au+Au y^s^N — 200 
AGeV and set to 42 mb Ta is the nuclear thickness 
function of nucleus A, 



Ta{s) 



dzpA{z, s). 



(11) 



where pa{z, s) is given by a Woods-Saxon parameteriza- 
tion of nuclear density, 

1 



PA{r) = po 



1 + Q{r~RA)/a ■ 



(12) 



In Eq. (|12|l parameters a, Ra, Po are 0.54 fm, 6.38 and 
0.1688, respectively On the other hand, the distri- 
bution of the number of binary collisions is given by 



^Nbc 



<yTA{bA)TB{bB). 



(13) 



Then total W{x,y;b) = w^!-f^ + (1 - w) ^ where 
w is the weight factor for binary scaling and is set to 
0.6, again utilizing a comparison of experimental data of 
single particle spectra to our model calculations. Fig- 
ures 121 and 13 show initial energy density in the longitu- 
dinal direction and on the transverse plane for Au-f Au 
y/SNN = 200 GeV central collisions for the case of a hy- 
brid hydro-|-micro calculation. 
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FIG. 2: Initial energy density in longitudinal direction at 



Au+Au ^/sNN ~ 200 GeV central collisions {h—2.4 fm) 
in the case of Hydro + UrQMD. 

b=2.4 fm - 




FIG. 3: Initial energy density on the transverse plane at 
Au+Au ^/sNN ~ 200 GeV central collisions {h—2.4 fm) 
in the case of Hydro + UrQMD. 

As a starting point we set initial longitudinal flow to 
Bjorken's scaling solution Q and neglect initial trans- 
verse flow. This is the simplest initial flow profile which 
will serve as basis for further investigation. For example, 
Kolb and Rapp discussed the possibility of existence of 
initial transverse flow which improves the results for Pt 
spectra and reduces the anisotropy in flow however, 
at the expense of introducing an additional parameter. 
Utilizing a parameterized evolution model, it has been 
pointed out that a Landau-type initial longitudinal com- 
pression and re-expansion of matter is favorable for the 
description of Hanbury-Brown Twiss (HBT) correlation 
radii ;5S]. This suggests that HBT analyses may be a 
sensitive tool for the determination of the initial longitu- 
dinal flow distribution. 

In Tab. Illll parameters with which we reproduce single 
particle spectra at RHIC are listed. The parameters 770 



and (Tri in longitudinal direction are determined mainly 
by the hadron rapidity distributions and do not strongly 
afl[ect the transverse momentum distributions. The pa- 
rameters for the initial conditions need to be optimized 
separately for the pure hydro and hydro + UrQMD cal- 
culations. A comparison between the two sets of initial 
conditions and possible physics implications can be found 
in Sec. IIIII Comparing the initial energy density of our 
purely hydrodynamic calculation to Ref. 57] in which 
an initial condition was parameterized in terms of the 
entropy- and baryon number density (weighted according 
to the wounded nuclear model by 75% and to the binary 
collision model by 25 %) we find that our calculation 
requires a significantly higher initial energy- or entropy 
density: the maximum value of the entropy density at Au 
-I- Au yJ/vTv = 200 GeV in central colHsion in [53 is 110 
fm~'^, whereas our value of the maximum entropy den- 
sity is 176 fm^^. The difference between the two values 
is most likely due to the treatment of resonance decays 
- these are explicitly treated in but are neglected 
in our case. Note that we consider our purely hydro- 
dynamic calculation solely as a baseline against which 
we can compare the full hydro+micro model, which of 
course contains a proper treatment of resonance decays. 
A third hydrodynamic model implementation described 
in uses an initial condition purely based on binary 
collisions. The maximum value of energy density in that 
calculation (which also treats resonance decays) is found 
to be ~ 45 GeV/fm'^, similar to the value we obtain for 
our hydro+micro calculation. 





TO 

(fm/c) 


£max 

(GeV/fm^) 


(fm-«) 


Vo 


(7,, 


pure hydro 


0.6 


55 


0.15 


0.5 


1.5 


hydro + UrQMD 


0.6 


40 


0.15 


0.5 


1.5 



TABLE III; Parameters for initial conditions of pure hydro 
and hybrid model. 



D. Hadronization and the transition to microscopic 
dynamics 

Having specified the initial conditions on the r = Tj 
hypersurface and the EoS, the hydrodynamical solution 
in the forward light-cone is determined uniquely. We 
assume that a freezeout process happens when a temper- 
ature in a volume element of fluid is equal to a freezeout 
temperature Tf in the pure 3-D hydrodynamic model. In 
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a hybrid model the transition from macroscopic to micro- 
scopic dynamics takes place at a switching temperature 
Tsw- The feezeout and switching temperatures respec- 
tively can be treated as parameters and determined by 
comparison to experimental data on single particle spec- 
tra. 

Due to our use of Lagrangian hydrodynamics, grid 
points move along the flux of the fluid and are not rep- 
resented by a fixed coordinate-space lattice any more. 
Therefore it is non-trivial to estimate of the number 
of particles flowing through the freezeout hypersurface. 
Here we start from a simple case [s^: suppose that the 
number of particles N{t) exists in the enclosed volume 
that is bounded by a closed surface S{t) at time r. N{t) 
is given by 



N{t) 



cfrn{r, r), 



(14) 



n(r) 



where n(r , r) is particle number density. At time t -|- 6t, 
the number of particle N{t) changes to 



dN 



d-^r 



dn 



1 



f2(r) 



dr dr 



d^rn{r, r), 



(15) 



sn 



where the volume varies to fl + SQ. Utilizing current 
conservation, {j is a current of particle), Ea. H15|) 

is rewritten as 

^ = - / d^sj ■n + [ d^s^n, (16) 
dr Js{t) Js{t) dr 

where n is the normal vector of surface element d^s and 
d( is the distance between the surface of ^(t) and that 
of fl{T) + Sfl. In Eq. dN/dr is the number of parti- 
cles which cross the surface S{t) during dr. Then total 
number of particles through the hypersurface S which is 
the set of surfaces {S'(r)}, 



N = / fda^, 



(17) 



where ~ n, dao = d^r, dcr = dr ■ d'^sn. If we write 

.„ gh 1 



E (27r)3exp[(P^w''-/i/)/r/]±l 



(18) 



for the current in Eq. H17|l . we obtain the Cooper- Frye 
formula [g^l 



E 



dN 



E 



9h 



(19) 



where gh is a degeneracy factor of hadrons and Tf and 
/i f are the freezeout temperature and chemical potential. 



In other words we obtain dcr^ by estimating the normal 
vector on the freezeout hypersurface S. Using Eq. H19|l . 
we then can calculate all particle distributions. 

In order to pass on the distribution H19|l in the macro- 
scopic model to the microscopic model, we first calculate 
the multiplicities Ni for each particle species i, by inte- 
gration of Eq. (jl9|l in which Tf and /j,/ are changed to 
the switching temperature Tsw and chemical potential 
Msw over space-time (t, r, 77) and momentum space (Pt, 
y). Ni is rounded to an integer value since the hadronic 
transport model described in the next section deals with 
real particle degrees of freedom. The distribution (|19(l 
divided by Ni is used as probability distribution to ran- 
domly generate the space-time and momentum-space co- 
ordinates for Ni hadrons of species i. This distribution 
serves as an input for the hadronic transport model for a 
single event. The sampling procedure can be repeated to 
generate a sequence of events as starting points for the 
microscopic calculation. Each event-sampling produces a 
different set of space-time and momentum coordinates as 
input of the microscopic model, however, the total mul- 
tiplicity for each species i remains constant for a given 
hydrodynamic switching hyper-surface. In a more real- 
istic calculation, particle number fluctuations should be 
taken into account as well. The reverse process, absorb- 
ing microscopic particles into the hydrodynamic medium 
is being neglected - it has been shown in that for 
rapidly expanding systems these contributions are negli- 
gible. 



E. Microscopic dynamics: the UrQMD approach 

The ensemble of hadrons generated accordingly is then 
used as initial condition for the microscopic transport 
model Ultra-relativistic Quantum Molecular Dynamics 
(UrQMD) [li|. The UrQMD approach is closely related 
to hadronic cascade |^ , Vlasov-Uehling-Uhlenbeck [b^ 
and (R)QMD transport models 0|- We shall describe 
here only the part of the model that is important for 
the application at hand, namely the evolution of an ex- 
panding hadron gas in local equilibrium at a temperature 
of about Tc ~ 160 MeV. The treatment of high-energy 
hadron- hadron scatterings, as it occurs in the initial stage 
of ultrarelativistic collisions, is not discussed here. A 
complete description of the model and detailed compar- 
isons to experimental data can be found in [l^ . 

The basic degrees of freedom are hadrons modeled as 
Gaussian wave-packets, and strings, which are used to 
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model the fragmentation of high-mass hadronic states via 
the Lund scheme . The system evolves as a sequence 
of binary collisions or 2 — iV-body decays of mesons, 
baryons, and strings. 

The real part of the nucleon optical potential, i.e. a 
mean-field, can in principle be included in UrQMD for 
the dynamics of baryons (using a Skyrme-type interac- 
tion with a hard equation of state). However, currently 
no mean field for mesons (the most abundant hadrons in 
our investigation) are implemented. Therefore, we have 
not accounted for mean-fields in the equation of motion of 
the hadrons. To remain consistent, mean fields were also 
not taken into account in the EoS on the fluid-dynamical 
side. Otherwise, pressure equality (at given energy and 
baryon density) would be destroyed. We do not expect 
large modifications of the results presented here due to 
the effects of mean fields, since the "fiuid" is not very 
dense after hadronization and current experiments at 
SIS and AGS only point to strong medium-dependent 
properties of mesons (kaons in particular) for relati vely 
low incident beam energies {Eiab < 4 GeV/nucleon) [6^. 
Nevertheless, mean fields will have to be included in the 
future. 

Binary collisions are performed in a point-particle 
sense: Two particles collide if their minimum distance 
d, i.e. the minimum relative distance of the centroids of 
the Gaussians during their motion, in their CM frame 
fulfills the requirement: 



d < do 



CTtot 



0-tot = cr(^/s, type) 



(20) 



The cross section is assumed to be the free cross section 
of the regarded collision type {N --N,N~A, tt~N...). 

The UrQMD collision term contains 53 different 
baryon species (including nucleon, delta and hyperon 
resonances with masses up to 2 GeV) and 24 differ- 
ent meson species (including strange meson resonances) , 
which are supplemented by their corresponding anti- 
particle and all isospin-projected states. The baryons and 
baryon-resonances which can be populated in UrQMD 
are listed in tabled the respective mesons in table ITTl - 
full baryon/antibaryon symmetry is included (not shown 
in the table) , both, with respect to the included hadronic 
states, as well as with respect to the reaction cross sec- 
tions. All hadronic states can be produced in string de- 
cays, s-channel collisions or resonance decays. 

Tabulated and parameterized experimental cross sec- 
tions are used when available. Resonance absorption, 
decays and scattering are handled via the principle of 



detailed balance. If no experimental information is avail- 
able, the cross section is either calculated via an One- 
Boson-Exchange (OBE) model or via a modified additive 
quark model which takes basic phase space properties 
into account. 

In the baryon-baryon sector, the total and elastic 
proton-proton and proton-neutron cross sections are well 
known 67] . Since their functional dependence on ^ysNN 
shows a complicated shape at low energies, UrQMD uses 
a table-lookup for those cross sections. However, many 
cross sections involving strange baryons and/or reso- 
nances are not well known or even experimentally acces- 
sible - for these cross sections the additive quark model 
is widely used. 

As we shall see later, the most important reac- 
tion channels in our investigation are meson-meson and 
meson-baryon elastic scattering and resonance formation. 
For example, the total meson-baryon cross section for 
non-strange particles is given by 



hot Ws) 



R=A,N 



{jB,mBjM, TUmW Jr, Mr) 



^Sr + I 



' {2Sb + 1){2Sm + 1) 

'PCMS (Mk- V^)2 + % 



(21) 



with the total and partial y^-dependent decay widths 
Ttot and Tr^mb- The full decay width Ttot{M) of a 
resonance is defined as the sum of all partial decay widths 
and depends on the mass of the excited resonance: 



TtotiM) 



E 

br={i,j} 



r.,,(M) 



(22) 



The partial decay widths ri,j{M) for the decay into the 
final state with particles i and j is given by 

21 + 1 



r.,,(M) 



Mr f {p^AM)) 



M 



1.2 



H-0.2 



\{P^AMR))) 



(23) 



here Mr denotes the pole mass of the resonance, T^j^ its 
partial decay width into the channel i and j at the pole 
and I the decay angular momentum of the final state. 
All pole masses and partial decay widths at the pole 
are taken from the Review of Particle Properties 67]. 
Fi J- (M) is constructed in such a way that Vij {Mr) — V^^ 
is fulfilled at the pole. In many cases only crude estimates 
for F^^ are given in jg^ - the partial decay widths must 
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then be fixed by studying exclusive particle production 
in elementary proton-proton and pion-proton reactions. 
Therefore, e.g., the total pion-nucleon cross section de- 
pends on the pole masses, widths and branching ratios 
of all N* and A* resonances listed in tabled Resonant 
meson-meson scattering (e.g. 7t + tt—^pottt + K—>- K*) 
is treated in the same formalism. 

In order to correctly treat equilibrated matter (we 
repeat that the hadronic matter with which UrQMD is 
being initialized in our approach is in local chemical and 
thermal equilibrium), the principle of detailed balance is 
of great importance. Detailed balance is based on time- 
reversal invariance of the matrix element of the reaction. 
It is most commonly found in textbooks in the form: 



pi 



(24) 



with g denoting the spin-isospin degeneracy factors. 
UrQMD applies the general principle of detailed balance 
to the following two process classes: 

1. Resonant meson-meson and meson-baryon interac- 
tions: Each resonance created via a meson-baryon 
or a meson-meson annihilation may again decay 
into the two hadron species which originally formed 
it. This symmetry is only violated in the case of 
three- or four-body decays and string fragmenta- 
tions, since N-body collisions with (N> 2) are not 
implemented in UrQMD. 

2. Resonance- nucleon or resonance- resonance inter- 
actions: the excitation of baryon-resonances in 
UrQMD is handled via parameterized cross sections 
which have been fitted to data. The reverse reac- 
tions usually have not been measured - here the 
principle of detailed balance is applied. Inelastic 
baryon-resonance de-excitation is the only method 
in UrQMD to absorb mesons (which are bound in 
the resonance). Therefore the application of the 
detailed balance principle is of crucial importance 
for heavy nucleus-nucleus collisions. 

Equation l|24|) . however, is only valid in the case of sta- 
ble particles with well-defined masses. Since in UrQMD 
detailed balance is applied to reactions involving reso- 
nances with finite lifetimes and broad mass distributions, 
equation H24() has to be modified accordingly. For the 
case of one incoming resonance the respective modified 
detailed balance relation has been derived in 'gs']. Here, 
we generalize this expression for up to two resonances in 
both, the incoming and the outgoing channels. 



The differential cross section for the reaction (1,2)^ 
(3 , 4) is given by: 



4=3 



here the pi in the (5-function denote four-momenta. The 
(5-function ensures that the particles are on mass-shell, 
i.e. their masses are well-defined. If the particle, how- 
ever, has a broad mass distribution, then the (5-function 
must be substituted by the respective mass distribution 
(including an integration over the mass): 

|■^|^^d^^^ 



daft 



dm 



P34- 



^=3 {m~M^f +Ty4 2n 



Incorporating these modifications into equation (|24|l and 
neglecting a possible mass-dependence of the matrix ele- 
ment we obtain: 

(pI^) (25i + 1)(25i + 1) 



(26) 



da|| 

dn 



(pI,) (2^3 + 1) (2^4 + 1) 

■■'+ J^34 

(ji"^ij2m2||JM)— y 

J=J- 



(27) 



Here, Si indicates the spin of particle i and the summa- 
tion of the Clebsch-Gordan-coefficients is over the isospin 
of the outgoing channel only. For the incoming channel, 
isospin is treated explicitly. The summation limits are 
given by: 

J_ = max(|ji-j2|,b3-j4|) (28) 
J+ = min(ji + j2,j3 + j4) • (29) 

The integration over the mass distributions of the reso- 
nances in equation 127|) has been denoted by the brackets 
0, e.g. 

pIa ^ (pIa) 

PcMsiVs^ TO4) ^3(7713) ^4(7714) dTO3 drn4 

with the mass distribution Ar{m) given by a free Breit- 
Wigner distribution with a mass-dependent width ac- 
cording to equation (|22|l : 

1 r(TO) 



Ar{'m) 



N {rrir - m)2 -hr(m)2/4 



(30) 



with 



lim Ar{m) = 6{mr — m) , 

and the normalization constant 

00 

r(m) 



N 



{rrir — m)2 -|- r(m)2/4 



■ dm . 



(31) 
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Alternatively one can also choose a Breit-Wigner distri- 
bution with a fixed width, the normalization constant 
then has the value N = 2-k. 

The most frequent applications of equation H27|l 
in UrQMD are the processes A1232 N N N and 
A1232 A1232 N N. 



III. RESULTS 

A. Modeling of the entire reaction dynamics in 
three-dimensional hydrodynamics 

The purpose of this section is to establish a baseline 
against which we can later compare the results of our hy- 
brid macro-f micro model. In addition we wish to demon- 
strate that our novel implementation of the 3+1 dimen- 
sional relativistic hydrodynamic model, utilizing a La- 
grangian grid, is capable of reproducing the results of 
previous hydrodynamic calculations. Many of these pre- 
vious calculations found in the literature have focused on 
one or two particular observables - here we wish to con- 
duct a consistent analysis of all relevant single-particle 
distributions which can be addressed by a hydrodynamic 
calculation. 

In order to obtain qualitative results which are com- 
parable to experimental data, we have to determine the 
parameters of the initial conditions. This is accom- 
plished by adjusting the parameters to fit the single par- 
ticle spectra for the most central Au+Au collisions at 
y/sWN = 200 GeV. 

First we show the behavior of isentropic trajectories 
in the T — fi plane for Au+Au ^/Tmn = 200 GeV cen- 
tral collisions in Fig. 21 The dotted line stands for 
the phase boundary between the QGP and the hadronic 
phase (Note that due to small baryochemical potentials, 
the phase boundary is an almost flat line at Tc — 160 
MeV). Apart from the central cell, we also investigate 
the isentropic trajectory of a cell close to the surface of 
the initially produced QGP. Whereas the isentropic tra- 
jectory of the central cell located at (0, 0, 0) starts in the 
QGP phase (solid line), the cell at the initial surface of 
the QGP (dashed line) only exhibits an evolution from 
the mixed phase to the hadronic phase. Both trajectories 
are terminated at freezeout temperature, Tf = 110 MeV. 

Figure El shows the Pt spectra of tt, K and p in Au + 
Au at y/sNN = 200 GeV for central collisions. Our cal- 
culation succeeds in reproducing the n spectra measured 
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FIG. 4: The isentropic trajectories on T ~ jj, plane. Solid 
(Dashed) line stands for the behavior of sell which is located 
at {x,y,rj) = (0,0,0)(= (0,6,0)) at initial time. Dotted line 
represents phase boundary. 



by the PHENIX collaboration '69"! up to Pt 2 GeV. 
However, due to the model assumption of chemical equi- 
librium up to the (low) kinetic freezeout temperature, we 
fail to obtain the correct normalization and hadron ra- 
tios, even though the shape of the Pt spectra of p and 
multistrange baryons (shown in Fig. (BJ is close to experi- 
mental data. In order to obtain the proper normalization 
for the p spectra and match the experimentally measured 
hadron ratios we adopt a procedure outlined in Ref . |7nj | , 
which renormalizes the Pt spectra using the p to tt ra- 
tio at the critical temperature to fix the normalization of 
the proton spectra. It is straightforward to extend this 
procedure to hyperons and multi-strange baryons as well, 
even though we choose to show the real, unrenormalized, 
result for the multi-strange baryons in Fig. Elto elucidate 
the situation prior to renormalization. 

The need for renormalizing the p spectra suggests that 
the assumption of a continuous chemical equilibrium un- 
til kinetic freezeout is not realistic and that an improved 
treatment of the freeze-out process is required. One 
method to deal with the separation of chemical and ther- 
mal freeze-out is the partial chemical equilibrium model 
(PCE) [fill I57I FH: below a chemical freeze-out temper- 
ature Tch one introduces a chemical potential for each 
hadron whose yield is supposed to be frozen-out at that 
temperature. The PCE approach can account for the 
proper normalization of the spectra, however, it fails to 
reproduce the transverse momentum and mass depen- 
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dence of the elliptic flow |72|- In section UlI Bl we shall 
utilize our hybrid hydro+micro model to decouple chem- 
ical and kinetic freeze-out. In these hybrid approaches 
[isl Q, ^3 the freeze-out occurs sequentially as a result 
of the microscopic evolution and flavor degrees of freedom 
are treated explicitly through the hadronic cross sections 
of the microscopic transport, leading to the proper nor- 
malization of all hadron spectra. 
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FIG. 5: Pt spectra for -n^ , , p at central collisions with 
PHENIX data [g^ . For proton we rescale our result using the 
ration at chemical temperature (See the text). 
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FIG. 6: Pt spectra for multistrange baryons at central colli- 
sions with STAR data fr^ - In this calculation (pure hydro) 
the additional procedure for normalization is not done. 



In Fig. 13 the centrality dependence of Pt spectra for 
TT"*" is shown. The impact parameters are set to b =2.4, 
4.5, 6.3, 7.9 fm corresponding to 0-6, 10-15, 15-20, 20- 



30 % centrality, respectively. These values are estimated 
via the number of nucleon-nucleon binary collisions and 
the number of participant nucleons in Ref. [s^. The 
centrality dependence is determined simply by the col- 
lision geometry - no additional parameter is necessary 
for our finite impact parameter collision calculations in 
Sec. Ill (]l Our results are consistent with the experi- 
mental data 69] in the low transverse momentum region 
(0 < Pt < 1 GeV) for all centralities. We observe that in 
peripheral collisions the difference between experimental 
data and our calculations appears at lower Pt compared 
to central collisions. This difference is indicative of the 
diminished importance of the soft, collective, physics de- 
scribed by hydrodynamics compared to the contribution 
of jet-physics in peripheral events. 
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FIG. 7: Centrality dependence of Pt spectra of 7r+ with 
PHENIX data 0. The Pt spectra at 10-15 %, 15-20 % 
and 20-30 % are divided by 5, 25 and 200, respectively. 



Figure IHl shows the centrality dependence of the pseu- 
dorapidity distribution of charged hadrons compared to 
PHOBOS data 73j. The parameters of our initial con- 
dition in the longitudinal direction (cr,j, 770) have been 
determined by fitting the data for the most central col- 
lisions, no additional parameter is needed to fix the ini- 
tial conditions for non-central collisions, since they re- 
sult solely from the collision geometry. The impact pa- 
rameters for our calculation are set to b =2.4, 4.5, 6.3 
and 7.9 fm corresponding to 0-6 %^-15 %, 15-25 % 
and 25-35 % centrality, respectively [73- We find good 
agreement of our results with the experimental data not 
only at mid rapidity but also at forward and backward 
rapidities, suggesting that the pure 3-D hydrodynamic 
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model can explain the charged hadron multiplicity dis- 
tribution in a wide range of rapidity. This observa- 
tion is consistent with results of other pure hydrody- 
namic models |2ol Isil utilizing different computa- 
tional methods. However, as we show it later, the hydro- 
dynamic calculation overestimates the elliptic flow at for- 
ward/backward rapidity, indicating that the rapidity dis- 
tribution is rather insensitive to the details of the expan- 
sion dynamics. Having fixed all parameters of our initial 
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FIG. 8: Centrality dependence of pseudorapidity distribution 
with PHOBOS data 1731 . Impact parameters in the hydrody- 
namic model are 2.4, 4.5, 6.3, 7.5 fm from central to peripheral 
collisions. 

condition utilizing the Pt spectra and (pseudo)rapidity 
distributions in the most central collisions we now apply 
our hydrodynamic model to other physical observables 
using these parameters. 

The elliptic flows as a function of Pt in 5-10 % and 
10-20 % most central collisions are shown in Fig. (a) 
and (b) together with STAR data [t^ at mid rapidity. 
We set the impact parameter in our hydrodynamic cal- 
culation to 4.5 fm (6.3 fm) for the 5-10 % (10-20 %) 
data. In the case of the 5-10 % most central events , 
we obtain reasonable agreement to the data for tt and 
p. On the other hand, in the 10-20 % centrality bin our 
hydrodynamical calculation overpredicts the elliptic flow 
compared to experimental data. Especially for protons, 
the deviation between calculation and experimental data 
is fairly large - this trend has already been observed in 

Q. 

Figure ^1 shows the elliptic flow as a function of r] in 
central (3-15 %) and mid central collisions (15-25 %). In 
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FIG. 9: Elliptic flow as a function of Pt with STAR data |7a| 
in centrality 5-10 % (a) and in centrality 10-20 % (b). Impact 
parameters in hydrodynamic calculations are 4.5 fm (a) and 
6.3 fm (b). 



both cases our hydrodynamic model calculations over- 
estimate the elliptic flow at forward and backward ra- 
pidities, similar to results shown in Ref. |51i|. At large 
forward and backward rapidities the assumptions of a 
perfect hydrodynamic model such as local equilibrium, 
vanishing mean free path and neglect of viscosity effect 
apparently are no longer valid. The deviations at forward 
and backward rapidities between experimental data and 
calculated results increase with the impact parameter, in- 
dicating a decrease of the volume in which hydrodynamic 
limit is achieved. 

Summarizing this section, we have applied our ideal 3D 
RFD model to Au+Au collisions at the top RHIC energy. 
A set of parameters for the initial conditions has been 
determined which allows for the simultaneous description 
of TT, K and p Pt spectra, the charged hadron rapidity 
distribution and the Pt as well as rapidity dependence 
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FIG. 10: Elli ptic flow as a function of r/ with PHOBOS experi- 
mental data [Tg for central (3-15 %) and mid central collisions 
(15-25 %). Impact parameters are set to 4.5 (central) and 6.3 
(mid central) fm, respectively. 



of the elliptic flow coefficient V2 for tt and p. Without 
any additional parameters the 3D RFD model is capable 
of describing the centrality dependence of the Px spectra 
and charged hadron rapidity distribution as well. 

However, our analysis also has demonstrated a couple 
of deficits of the ideal 3D RFD approach - many of which 
are already well-known and have been discussed in the 
literature before [l5[l3| : 

• the centrality dependence of the elliptic flow coef- 
ficient V2 as a function of Pt is not well described, 

• the width of the V2 vs. rj distribution is too broad 
compared to data, 

• p and multistrange particle spectra need to be nor- 
malized by hand in order to account for the sepa- 
ration of chemical and kinetic freeze-out, 

• the hydrodynamic approach is only of limited ap- 
plicability for small systems (i.e. large impact pa- 
rameters) and large Pt (hard physics). 

In the following section, we shall use the 3D RFD cal- 
culation as a baseline to determine the effects of an im- 
proved treatment of the hadronic phase in the framework 
of the hydro-|-micro approach. 



B. Application of the hydro+micro approach 

As in the previous subsection, we first determine the 
parameters of the initial conditions in Sec. lIICJl bv fitting 
the single particle spectra in the most central centrality 
bin. Table ITTll shows the parameters for both, the pure 
3D RFD initial condition as well as the hydro-f micro ini- 
tial condition. The main difference we find between the 
two is in the maximum value of initial energy density. 
The large difference between the two initial energy den- 
sities can be explained by our omission of resonance de- 
cays in the pure hydrodynamic calculation and to a lesser 
extent by the dissipative corrections present in the hy- 
dro+micro calculation. Both effects result in additional 
particle production, leading to a smaller initial energy- 
and entropy-density necessary to describe the final par- 
ticle multiplicities. 

We set the switching temperature Tgw to 150 MeV. 
This allows for a brief period of hydrodynamic evolu- 
tion in the hadronic phase to account for multi-particle 
collisions which can occur at large densities and temper- 
atures in the hadronic phase close to the phase-boundary 
The dependence of hadronic observables on Tswj 
such as the mean transverse momentum {Pt) of the dif- 
ferent hadron species, has been investigated in [ifil l - our 
choice of 150 MeV conforms to the lower bound of the 
allowed range for Tsw- 

Figure ^2 shows the Pt spectra of 7r+ , K and p at 
y/SNN = 200 GeV central collisions. The most com- 
pelling feature, compared to the pure 3D RFD calcula- 
tion, is that the hydro-|-micro approach is capable of ac- 
counting for the proper normalization of the spectra for 
all hadron species without any additional correction as is 
performed in the pure hydrodynamic model. The intro- 
duction of a realistic freezeout process provides therefore 
a natural solution to the problem of separating chemical 
and kinetic freeze-out in a pure hydrodynamic approach. 

In Fig. I12lcentralitv dependence of Pt spectra of 7r+ is 
shown. The impact parameter for each centrality is de- 
termined in the same way as in the pure hydrodynamic 
calculation. The separation between model results and 
experiment appears at lower transverse momentum in pe- 
ripheral collisions compared to central collisions, just as 
in the pure hydrodynamic calculation. The 3D hydro -t- 
micro model does not provide any improvement for this 
behavior, since the hard physics high Pt contribution 
to the spectra occurs at early reaction times before the 
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FIG. 11: Pt spectra for tt , K and p at central collisions 
with PHENIX data 0|. The poins are NOT renormalized. 
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system has reached the QGP phase and is therefore nei- 
ther included in the pure 3D RFD calculation nor in the 
hydro-|- micro approach. 

Figure [TSl shows the centrality dependence of the pseu- 
dorapidity distribution of charged hadrons compared to 
PHOBOS data 0. Sohd circles stand for model re- 
sults and open circles denote data taken by the PHO- 
BOS collaboration [t^. The impact parameters are set 
toh = 2.4,4.5,6.3,7.9 fm for 0-6 %, 6-15 %, 15-25 % 
and 25-35 % centralities, respectively. Our results are 
consistent with experimental data over a wide pseudo- 



rapidity region. We observe a small deviation around 
1 77 1 3, which may be improved by tuning the param- 
eter cr,j (here we have chosen the same value as for the 
pure RFD calculation). There is no distinct difference 
between 3-D ideal RFD model and the hydro + UrQMD 
model in the centrality dependence of the psuedorapidity 
distribution, indicating that the shape of psuedorapidity 
distribution is insensitive to the detailed microscopic re- 
action dynamics of the hadronic final state. 
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FIG. 13: Centrality dependence of pseudorapi dity distribu- 
tion of charged particles with PHOBOS data |73ll . Impact 
parameters in a calculation are 2.4 (0-6 %), 4.5 (6-10 %), 6.3 
(10-15 %), 7.5 fm (25-35 %). 



In Fig. ^] we analyze the Pt spectra of multistrange 
particles. Our results show good agreement with experi- 
mental data for A, S, Vl for centralities 0-5 % and 10-20 
%. In this calculation the additional procedure for nor- 
malization is not needed. Recent experimental results 
suggest that at thermal freezeout multistrange baryons 
exhibit less transverse flow and a higher temperature 
closer to the chemical freezeout tem per ature compared 
to non- or single-strange baryons [t^IsoI . This behavior 
can be understood in terms of the flavor dependence of 
the hadronic cross section, which decreases with increas- 
ing strangeness content of the hadron. The reduced cross 
section of multi-strange baryons leads to a decoupling 
from the hadronic medium at an earlier stage of the re- 
action, allowing them to provide information on the prop- 
erties of the hadronizin g Q GP less distorted by hadronic 
final state interactions [sii Is^ l . It should be noted that 
the analogous behavior has already been observed in ex- 
periments at the CERN-SPS Isil. Later in this section 
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we will discuss the reaction dynamics of multi-strange 
baryons in greater detail by analyzing the baryon colli- 
sion number and freezeout time distributions as well as 
their collision rates. 
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FIG. 14: Pt spectra of multi-strange particles at centralities 
0-5 % and 10-20 % with STAR data 



In Fig. 1151 the mean transverse momentum (Pt) as a 
function of hadron mass is shown. Open symbols de- 
note the value at Tsw = 150 MeV, corrected for hadronic 
decays. Not surprisingly, in this case the (Pt) follow 
a straight line, suggesting a hydrodynamic expansion. 
However if hadronic rescattering is taken into account 
(solid circles) the (Pt) do not follow the straight line any 
more: the (Pr) of pions is actually reduced by hadronic 
rescattering (they act as a heat-bath in the collective 
expansion), whereas protons actually pick up additional 
transverse momentum in the hadronic phase. RHIC data 
by the STAR collaboration is shown via the solid trian- 
gles - overall the proper treatment of hadronic final state 
interactions significantly improves the agreement of the 
model calculation with the data. 

Let us now investigate the effect of resonance decays 
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and hadronic rescattering on the pion and baryon trans- 
verse momentum spectra: Figure [TBI shows the Pt spec- 
trum for 7r+ at Tgw = 150 MeV (solid line, uncorrected 
for resonance decays) as well as the final spectrum after 
hadronic rescattering and resonance decays, labeled as 
Hydro+ UrQMD (solid symbols). In addition the open 
symbols denote a calculation with the resonance de- 
cay correction performed at Tgw, which we label as Hy- 
dro + decay. The difference between the solid line and 
open symbols therefore directly quantifies the effect of 
resonance decays on the spectrum, which is most dom- 
inant in the low transverse momentum region Pt < 1 
GeV Furthermore, the comparison between open sym- 
bols and solid symbols quantifies the effect of hadronic 
rescattering: pions with Pt > 1 GeV lose momentum via 
these final state interactions, resulting in a steeper slope. 

Figure [T7I shows a likewise analysis for baryons (p. A, 
S and ri). We note that in contrast to SPS energies |15l |. 
there is very little effect on the spectra due to hadronic 
rescattering - even for protons which rescatter 8-10 
times at mid-rapidity. Only at low transverse momenta 
the multiple scatterings of the protons (predominantly 
with pions) manifests itself in a slight fiattening of the Pt 
distribution of the protons, giving rise to a slight increase 
in their radial flow. This phenomenon has been discussed 
in Ref. 113] and is commonly referred to as 'pion wind'. 

Figure IT51 shows the distribution of the number of col- 
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FIG. 17: Pt spectra of baryons from hydro-|-decay (open sym- 
bols) and hydro -|- UrQMD (solid symbols) at centrality 0-5 
%. 



lisions which particles suffer in the hadronic phase at 
b = 2.4, 4.5, 6.3 fm. The distributions for protons and 
lambdas are very broad, indicating a large amount of 
rescattering taking place in the hadronic phase. Around 
mid-rapidity protons rescatter on average 10 times for 
central and semi-central impact parameters and lambdas 
rescatter 7-8 times (See also Fig. ^]). However, the av- 
erage number of collision of multistrange baryons is less 
than half of those for protons and A's, as is to be ex- 
pected due to the decrease of the hadronic cross section 
with increasing strangeness content of the hadron. We 



note that the collision number distributions change sig- 
nificantly as a function of centrality - for large impact 
parameters the high-ricoii tail of the distribution exhibits 
a much steeper drop as a function of ricou- 
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FIG. 18: The impact parameter dependence of distribution 
of number of collisions for p, A, H, Q. 

The pseudo-rapidity dependence of the number of 
hadronic rescatterings for different baryon species is an- 
alyzed in Fig. 1191 which shows the number of collisions 
of p, A, S and as a function of j] at b = 2.4, 4.5 and 
6.3 fm. The distributions appear to be similar to that of 
the particle yield psuedorapidity distribution. At midra- 
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pidity we find a plateau region extending from r] = —3 to 
3, followed by a steep drop-off to forward and backward 
rapidities. The flavor dependence of the average collision 
numbers is again clearly seen, even though we would like 
to point out that the shapes of the different distributions 
is very similar. The large plateau region indicates the ra- 
pidity domain in which interacting matter can be found 
and in which the application of thermodynamic concepts 
is viable. 
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Figure QUI depicts the freezeout time distribution of tt, 



K, p, A, S and at midrapidity for b = 2.4, 4.5, 6.3 fm. 
The freezeout time distribution of tt and K has a peak 
around t ^ 19 {< 20 fm), while the peak of p and A 
freezeout time distribution is shifted to later times. We 
also note that the freezeout time of multistrange baryons 
is much smaller than that of other particles. The freeze- 
out times are in general determined by the amount of 
hadronic rescattering suffered by the different hadron 
species - our analysis indicates that (a) hadronic freeze- 
out is strongly species dependent, and (b) even for a par- 
ticular species, the freeze-out distribution is broad and it 
is therefore nearly impossible (or at least extremely am- 
biguous) to define a precise freeze-out time for a given 
hadron species. 

This figure again supports the finding that multi- 
strange baryons are produced with few interactions right 
after hadronization and freeze-out early 0, 0| , indica- 
tions of which have been observed by the STAR collabo- 
ration [t^ EBI ■ Given the broad freeze-out time distribu- 
tions it is very difficult to quantify the overall duration of 
the heavy-ion reaction - a possible criterion would be the 
drop-off of the number of freezing out particles per unit 
time and rapidity below 1 - which would put the overall 
duration of the reaction to approximately 30 fm/c. 

In Fig.l^we study the effect hadronic rescattering has 
on the duration of the freeze-out process by comparing 
a calculation terminated at Tgw without hadronic rescat- 
tering (open symbols) to one including the full hadronic 
final state interactions (solid symbols). If we terminate 
at Tgw, most hadrons freezeout around 10 fm/c, reflect- 
ing the lifetime of the dcconfincd phase in our calculation 
(the tails of the distribution stem from the decays of long- 
lived resonances). The inclusion of hadronic rescatter- 
ing shifts the peak of the freezeout distribution to larger 
freezeout times (tj ^ 20 — 30), providing us with an esti- 
mate on the lifetime of the hadronic phase around 10-20 
fm/c. Note that this estimate is subject to the same sys- 
tematic uncertainties discussed previously in the context 
of the overall lifetime of the system. 

The findings discussed in the context of the previous 
figure for pions and kaons are confirmed by analyzing 
baryons in the same fashion, which is shown in Fig. 1221 
here the top frame contains the analysis terminated at 
Tsw and the bottom frame contains the calculation in- 
cluding full hadronic rescattering. 

Figure 1231 shows collision rates for meson-meson, 
meson-baryon and baryon-baryon collisions at mid- 
rapidity and an impact parameter of b = 2.4 fm. We 
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r2 at mid rapidity for b = 2.4, 4.5, 6.3 fm. 



find that the reaction dynamics in the hadronic phase 
is dominated by meson-meson and to a lesser degree by 
meson-baryon interactions. Baryons essentiaUy propa- 
gate in a medium dominated by mesons - a situation 
very difl^erent from the regime at AGS and lower SPS 
energies. 

Finally, let us move to the analysis of elliptic flow in the 
hydro-|-micro approach: First we show the elliptic flow 
coefiicient V2 as a function of Pt for 7r+ at centralities 
5-10 % and 10-20 % in Fig. El Open symbols stand 
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FIG. 21: The freezeout time distribution of mesons for hy- 
dro+decay (open symbols) and hydro -f UrQMD (solid sym- 
bols) at mid rapidity in the case of central collision. 



for experimental data and solid symbols represent our 
calculations. We find that the hydro-|- micro approach 
is able to provide an improved agreement to the data 
compared to the purely hydrodynamic calculation (see 
Fig. for comparison). 

Figure [23 shows the elliptic flow coefficient V2 as a 
function of pseudo-rapidity rj for charged particles in cen- 
tral (3-15 %) (a) and mid central (15-25 %) (b) col- 
lisions. The solid line stands for a purely hydrody- 
namic calculation whereas the solid circles denote the 
hydro+micro approach and the solid triangles represent 
PHOBOS data 76]. We find that the dissipative effects 
contained in the hydro-|-micro approach significantly al- 
ter the shape of the V2 vs. rj curve, providing a far bet- 
ter agreement to the data - in particular for pseudo- 
rapidities away from mid-rapidity. The analysis shows 
that apparently dissipative effects increase towards pro- 
jectile and target rapidities, where the assumptions of 
ideal fluid dynamics break down earlier. Our analysis 
conflrms the findings of 0, who in addition also stud- 
ied the effect of a Color-Glass initial condition on the 
elliptic flow rapidity dependence. 

The question of how much elliptic flow develops dur- 
ing the deconfined phase vs. the hadronic phase is in- 
vestigated in Figs. 1211 and |23 It has been pointed out 
repeatedly, both in the framework of microscopic [s^ ^ 
well as as hydrodynamic analysis [3| that elliptic flow 
develops early on during the deconfined phase of the re- 
action and thus serves as a sensitive tool to the equa- 
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FIG. 22: The freezeout time distribution of baryons for hy- 
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tion of state of QCD matter. This picture has been 
confirmed by the experimentally observed quark number 
scaling V2 ^ nv2{l/nPT) of elliptic flow at intermediate 
transverse momenta Pt [s^ • 

In Fig. |2Slwe plot elliptic flow V2 as a function of Pt- 
The solid line stands for the pure hydro calculation, ter- 
minated at the switching temperature Tg^ and solid cir- 
cles denote the full hydro+micro calculation. Wc find 
that the QGP contribution to the elliptic flow depends 
on the transverse momentum - for low Pt nearly 100% of 
the elliptic flow is created in the QGP phase of the reac- 
tion, whereas the hadronic phase contribution increases 
to 25% at a Pt of 1 GeV/c. 

Figure shows the elliptic flow as a function of 77: 
the pure hydrodynamic calculation is shown by the solid 
curve, the hydrodynamic contribution at Ts„ is denoted 
by the dashed line and the full hydro-|- micro calculation is 
given by the solid circles, together with PHOBOS data 
(solid triangles). The shape of the elliptic flow in the 




FIG. 23: Collision rate for meson-meson (MM), meson-baryon 
(MB) and baryon-baryon (BB) at central collisions. The 
(dark) grey zone stands for qgp (mixed) phase, which is de- 
termined by the center cell. 
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FIG. 24: Elliptic flow as a function of Pt for tt^ at centralities 
5-10 % and 10-20 %. Open symbols stand for STAR data and 
solid symbols stands for our results. 



pure hydrodynamic calculation at Tgw is quite different 
from that of the full hydrodynamic one terminated at a 
freeze-out temperature of 110 MeV. Apparently the slight 
bump at forward and backward rapidities observed in the 
full hydrodynamic calculation develops first in the later 
hadronic phase, since it is not observed in the calcula- 
tion terminated at Tgw Evolving the hadronic phase in 
the hydro-|-micro approach will increase the elliptic flow 
at central rapidities, but not in the projectile and target 
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rapidity domains. As a result, the elliptic flow calcula- 
tion in the hydro+micro approach is closer to the exper- 
imental data when compared to the pure hydrodynamic 
calculation. 

IV. SUMMARY AND OUTLOOK 

In summary, we have introduced a hybrid macro- 
scopic/microscopic transport approach, combining a 
newly developed relativistic 3+1 dimensional hydrody- 
namic model for the early deconfined stage of the reaction 
and the hadronization process with a microscopic non- 



FIG. 27: Elliptic flow as a function of 77 of charged particles. 
Solid line and dashed lines stand for pure hydro calculation 
at the freezeout temperature (110 MeV) and the switching 
temperature, respectively. Experimental data by PHOBOS 
is shown with solid triangles. Open (Solid) circles stand for 
hydro -I- decay (hydro -I- UrQMD). 

equilibrium model for the later hadronic stage at which 
the hydrodynamic equilibrium assumptions are not valid 
anymore. Within this approach we have dynamically cal- 
culated the freezeout of the hadronic system, accounting 
for the collective flow on the hadronization hypersurface 
generated by the QGP expansion. We have compared 
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the results of our hybrid model and of a calculation uti- 
lizing our hydrodynamic model for the full evolution of 
the reaction to experimental data. This comparison has 
allowed us to quantify the strength of dissipative effects 
prevalent in the later hadronic phase of the reaction, 
which cannot be properly treated in the framework of 
ideal hydrodynamics. 

Overall, the improved treatment of the hadronic phase 
provides a far better agreement between transport cal- 
culation and data, in particular concerning the flavor 
dependence of radial flow observables and the collective 
behavior of matter at forward/backward rapidities. We 
find that the hadronic phase of the heavy-ion reaction 
at top RHIC energy is of significant duration (at least 
10 fm/c) and that hadronic freeze-out is a continuous 
process, strongly depending on hadron flavor and mo- 
menta. 

With this work we have established a base-line - both 
for the regular 3-1-1 dimensional hydrodynamic model as 
well as for the hybrid hydro-|-micro approach. In forth- 
coming publications we shall expand on this baseline by 
investigating the effects of a realistic lattice-QCD mo- 
tivated equation of state containing a tri-critical point 



and by performing an analysis of two particle correlations 
(HBT interferometry) . We also plan to use our model as 
the medium for the propagation of jets and heavy quarks 
and to study the modification of our medium due to the 
passage of these hard probes. 
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V. APPENDIX 

Equation in {T,x,y,ri) is written in the following 
explicit way. 
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+ 



1 ~ ~o 
1 ~ -2-2 



0, 



(32) 



where 7 = l/,/l - vl ~ - vl, Yl ^ l/21n((l + in the (r, a;, y, 77) coordinates is the addition of the fifth 
'y^)/(l - w^)). The diflFerence between the rclativistic hy- term in Eq. 
drodynamic equations in Cartesian coordinates and those 
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